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Abstract 

We present results of numerical simulations of the 3+1 dimensional Nambu - Jona-Lasinio (NJL) 
model with a non-zero baryon density enforced via the introduction of a chemical potential [i ^ 0. 
The triviality of the model with a number of dimensions d > 4 is dealt with by fitting low energy 
constants, calculated analytically in the large number of colors (Hartree) limit, to phenomenological 
values. Non-perturbative measurements of local order parameters for superfluidity and their related 
susceptibilities show that, in contrast with the 2+1 dimensional model, the ground-state at high 
chemical potential and low temperature is that of a traditional BCS superfluid. This conclusion is 
supported by the direct observation of a gap in the dispersion relation for 0.5 < ua < 0.85, which at 
fia = 0.8 is found to be roughly 15% the size of the vacuum fermion mass. We also present results 
of an initial investigation of the stability of the BCS phase against thermal fluctuations. Finally, 
we discuss the effect of splitting the Fermi surfaces of the pairing partners by the introduction of 
a non-zero isospin chemical potential. 
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I. INTRODUCTION 



The ground state of QCD with Nf > 3 flavors of light quark at low temperature T 
and asymptotically high density is thought to be a "color-flavor locked" (CFL) state in 
which SU(3) c (S>SU(iVy) L ®SU(iVj) R ®U(l)B symmetry is spontaneously broken by diquark 
condensation in the color anti-triplet channel to a diagonal SU(3)a, and which is thus 
simultaneously color superconducting, superfluid, and chirally broken {jH^. Since diquark 
condensation is thought to occur at the Fermi surface via a BCS mechanism, it is accurately 
described by a self-consistent gap equation in the asymptotic regime /x — ► oo, with /x the 
quark chemical potential or Fermi energy, where the QCD coupling g(/x) is weak 3(. However, 
the densities for which analytic predictions of weakly-coupled QCD can be trusted are much 
greater than the conditions, with /x ~ O(lGeV), which are likely to be physically realisable 
in our universe at the centres of compact stars j^]. 

In this regime we face the twin problems that QCD becomes strongly interacting, and 
that the most systematic way of computing its properties non-perturbatively, namely nu- 
merical simulation of lattice gauge theory, cannot help because of the "sign problem" ; the 
lack of positivity of the QCD Euclidean path integral measure with /x 7^ makes Monte 
Carlo sampling methods inoperable. Another complicating factor is that once the density 
becomes low enough to make the strange quark mass m s no longer negligible compared to 
/x, then other channels involving pairing between just u and d may be preferred, leading to 
a plethora of different possible ground states, including even crystalline examples BOB. 
Analytic approaches to these questions must either use some approximate non-perturbative 
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approach such as the instanton liquid [8(, or resort to phenomenological models of the strong 
interaction such as the Nambu - Jona-Lasinio (NJL) model 

Generic NJL models contain fermions, to be thought of either as quarks or baryons, 
self-interacting via contact four- Fermi terms 11]. In a Euclidean metric, the prototype is 
written in terms of isopinor fermion fields as 

Cnjl = </# + m)V - y [Wf ~ 0K75^) 2 ] , (1) 

which has SU(2)i®SU(2)^ chiral symmetry. In 3+1 dimensions the interaction strength g 2 
has mass dimension -2, implying that the model is non-renormalisable and must be defined 
in terms of some explicit ultra-violet scale A (see e.g. |l2|). Since the model has no gluons, 
it fails to reproduce the physics of confinement; however, for sufficiently large g 2 the model 



does exhibit spontaneous chiral symmetry breaking to SU(2)j, resulting in a physical or 
"constituent" fermion mass scale S ~ <? 2 A 3 ^> m. The phase structure of the model is most 
conveniently studied in the Hartree approximation, corresponding to retaining only those 
diagrams which, if the number of fermion degrees of freedom were multiplied by a factor iV 
and the coupling rescaled as g 2 /N, would remain at leading order in 1/JV. At low T it is found 
that for values of chemical potential fi c ~ E there is a transition, whose order depends on the 
details of the cutoff, in which chiral symmetry is restored and baryon density ub = (ipjoip) 
increases sharply from zero Q Q • For \i > \x c the N JL model thus describes a state 
resembling relativistic "quark matter." Since there is no short-range repulsion present, the 
stability of a phase in which both % and X are simultaneously non-zero, corresponding to 
"nuclear matter," is not a priori clear. 

The color superconducting solutions discussed in the first paragraphs have been obtained 
by solving the self-consistent "gap equation" for the smallest energy 2A required to excite 
a quasiparticle pair out of the ground state consisting of a filled Fermi sea. Solutions with 
A > 0, implying the instability of a sharp Fermi surface with respect to formation of a 
condensate of Cooper pairs, form the basis of the BCS description of superconductivity 
Such solutions are also characterised by a Cooper pair or diquark condensate (qq) ^ 
whose precise definition will be given below; here it suffices to identify it with the density 
of condensed pairs in the ground state. Since the NJL model is not a gauge theory the 
corresponding BCS ground state is not superconducting (analogous to a Higgs vacuum 
in particle physics), but rather a superfluid, in which the global U(1)b baryon number 
symmetry of is spontaneously broken by (qq) ^ 0, which is thus a true order parameter. 

Since to leading order in 1/N (qq) = A = 0, it is legitimate to query the approximate 
methods used to find such solutions. Fortunately in this case it is possible to employ nu- 
merical lattice methods, because as reviewed below the lattice-regularised NJL model has 
no sign problem. Initial studies have focused on the high density phase of the NJL model 
in 2+ld, in part for the obvious technical advantage of working in reduced dimensionality, 
and in part for the conceptual advantage that NJL2+1 has an ultra-violet stable fixed point 
and hence an interacting continuum limit, so that in principle we need not worry about 
the cutoff dependence of the results. Although the simulations reveal enhanced diquark 
pairing in the scalar isoscalar channel JitJ , there is no evidence for the expected BCS signal 
(qq) / 0, A ^ 0. Instead, the condensate vanishes non-analytically with external diquark 
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source strength j, and the fermion dispersion relation reveals a sharp Fermi surface and 
massless quasiparticles 18, 3]. The fitted values for the Fermi momentum and Fermi veloc- 
ity, however, do not match those of free field theory, and indeed are difficult to account for 
in orthodox Fermi liquid theory j^. The interpretation is that due to its low dimension- 
ality the system realises superfluidity in an unconventional Berezinskii-Kosterlitz-Thouless 
(BKT) mode first invoked to describe thin helium films 2l|; the massless modes due to 
phase fluctuations of the would-be condensate field remain strongly fluctuating, resulting in 
unbroken U(1)b symmetry but critical behaviour for all /i > /i c . Met ast ability of persistent 
flow is only revealed by the non-trivial response of the system to a symmetry breaking field 
with a twist imposed across the spatial boundary. 

For this reason in the current paper we present results of simulating NJL 3+1 , with the 
goal of exposing superfluid behaviour with the more conventional signal (qq) 7^ 0, A > 0. 
Our first results suggesting that (qq) 7^ for \i > fi c have appeared in |22|. The details 
of the lattice model, the simulation technique, and the observables monitored to expose 
this behaviour are reviewed below in Sec. El As already stated, in 3+ld the model is non- 
renormalisable and must be defined with a cutoff. The bare parameters g 2 and m must be 
chosen so that the model's properties match those of low energy QCD as closely as possible; 
in Sec. Ill CI we use a large- N expansion (here N is identified with the number of "colors" 
N c ) to calculate the vacuum fermion mass S and the constants f n and m^, and find that a 
reasonable matching can be made with an inverse lattice spacing a -1 ~ 700MeV. In Sec. IHII 
we present results for the phase structure of the model for low T, the diquark condensate 
and associated susceptibilities, and in Sec. HYI dispersion relations in the spin-| sector. For 
the first time there is evidence for a non-vanishing BCS gap from a lattice simulation, with 
A/Sq — 0.15, which translates to a gap of around 60MeV in physical units. This is consistent 
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2311 . In Sec. El we discuss finite volume 



with estimates from self-consistent approaches 
effects and finally in Sees. IVll and IV 111 present preliminary results for simulations both with 
T > and with small but non-zero isospin chemical potential /// 7^ 0, which has the effect 
of separating the u and d quark Fermi surfaces and hence discouraging ud condensation. 
Sec. IVIIII contains a brief summary. 
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II. LATTICE MODEL AND PARAMETER CHOICE 



A. Formulation and symmetries 

The lattice version of the NJL model we employ was first used in a study of chiral 



12| . The action, with lattice spacing a 



symmetry restoration at zero chemical potential in 
set to unity reads 

So = S ferm + S bos = Y,XxM[<5>] xyXy + ( x M*[$} xy ( y + 4 ^tr$t$ £ , (2) 

xy 9 x 

where Xi Xi C an d C are Grassmann- valued staggered isospinor fermion fields defined on the 
sites x of a hypercubic lattice, and $ = o + iix.f is a 2 x 2 matrix of bosonic auxiliary fields 
defined on the dual sites x. The kinetic operator M is given by 



M pq = \$ P q 
xy 2 



eM( Wo- e ^ S yx-0+ Yl vAx)(Syx+0-S y x-u) + 2m8 xy 

u=l,2,3 



+ ^S xy V {a{x)5™ + ie{x)fi(x).i**) (3) 
16 f-^ 

such that the parameters are bare fermion mass m, coupling g 2 (with mass dimension -2) 
and baryon chemical potential fi. The symbol (x, x) denotes the set of 16 dual sites adjacent 
to x. The Pauli matrices acting on isospin indices p, q are normalised so that trrj-r, = 28y. 
The phase factors r}„(x) = (— iyo+---+x^-i anc [ £ ^ _ ^_-iY 0+Xl+X2+Xi ensure that fermions 
with the correct Lorentz covariance properties emerge in the continuum limit, and that in 
the limit m — > the action © has a global SU(2)i®SU(2)/j invariance under 

X^(V e U + V V)x ; x^x(v e v^ + v u^) 
C (V e U* + V V*)C ; C ^ ((VeV^* + Votf*) 

$ y$C/t witn u,VeSV(2). (4) 

The even/odd projectors are given by V e / (x) = |(1 ± e(x)). In addition the action has a 
U(l)s invariance under baryon number rotations: 

X,C^e ia xX X,C^xXe~ ia . (5) 

The auxiliary fields a and 7f can be integrated out to yield an action written solely in 
terms of fermion fields x an d C m the continuum limit ad — > lattice fermion doubling 



leads to a physical content for the model of eight fermion species (four described by Xi f° ur 
by Q with self-interactions resembling those of the continuum NJL model (but with terms 
of the form xxCC differing in detail from those of the form XXXX) > resulting in an additional 
approximate U(4)®U(4) symmetry which is violated by terms of 0(a). Further details are 
given in [12| . In what follows we will somewhat arbitrarily refer to these degrees of freedom 
as colors and hence define N c = 8; this is in distinction to the explicit flavor or isospin 
symmetry which gives rise to Nf = 2. 

In order to examine issues associated with diquark pairing, it is necessary to introduce 
a source term into the action. A suitable term, invariant under (0J) but violating baryon 
number conservation (J5J), is: 

S source = £ \ (X^X, + £t 2 Q + \ (Xx^tl + , (6) 

x 

where in this paper we will consider j,j real and positive. With this addition the partition 
function follows from integrating over the Grassmann degrees of freedom: 

= J DxDxD(D(D®exp -{S ferm + S bos + S source ) 

= J D<$>Yi Nc l\A[<S>, 3 ,3})e- Sb °^ (7) 

with the antisymmetric matrix A(j,j), which acts on bispinors (x,X tr ) (and A*(— j, —j) on 
(CO), given by 

A= \ It • (8) 

\-M tr JT 2 J 

The reality of P£A = VdetA follows by notingthe identity t 2 Mt 2 = M* 0. In fact, PfA 
can also be shown to be positive definite jlfll l24j | , implying that Z can also be defined for 
a single power of the pfaffian, (i.e. with the fields (,( omitted), corresponding to iV c = 4. 
However, the smallest value of N c permitted by the exact hybrid Monte Carlo algorithm 



(see Sec. Ill Dl) is 8. 



B. Observables 



Here we list the various observables monitored in the simulations. For simplicity we 
restrict our attention to those constructed from the Xi X fields; their equivalents in the ( 
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sector are easily found. Angled brackets denote averages taken with respect to the measure 
defined by (JTj). 

The chiral condensate (xx) is given by 



, . 1 d\nZ 1 / / 1\ , , , 



To leading order in a large- iV c expansion the chiral condensate is proportional to the expec- 

2 

tation value of the scalar field S = (a) = \(xx)-> which in the chiral limit also defines the 
physical or constituent fermion mass. The baryon charge density ub per flavor is given by 

1 dlnZ ll ( e»6 yx+6 + e-»5 A \ 

-^^(H.v^v, r (10) 

In the diquark sector it is convenient to define operators 

qq±(x) = -jb&TzXx ± Xx^Xx] (ii) 

and corresponding source strengths j± — j ± J. The condensates are then 



1 dlnZ 1 / /±r 2 \ i i \ 
\ r 2 > 



V dj± AV 



In the continuum limit, the corresponding 
spin, 2 flavor and 4 color components are 



operators written in terms of spinors ?/>, ip with 4 



24] 



gg± oc V C 75 ® r 2 ® (C 75 )*^ ± ^C lh ® r 2 ® (C 75 )>* r , (13) 

with C the charge conjugation matrix satisfying C^^C^ 1 = — 7*. The condensate wavefunc- 
tion is thus scalar isoscalar, but with a non-trivial (and for us uninteresting) variation under 
"color" rotations, and manifestly antisymmetric as required by the exclusion principle. We 
will also consider the diquark susceptibilities 

x± = = 5>? ± (0)gg ± (*)) (14) 

which are constrained by identities analogous to the axial Ward identity for the pion prop- 
agator; a particularly useful form in the j_ — > limit reads 

X-b-- = ^. (15) 
3+ 



When U(l)s is spontaneously broken by the formation of a condensate (qq+) 7^ 0, therefore, 
gg_ interpolates the resulting Goldstone mode, whose masslessness as j + — > is guaranteed 
by (|15|) . Physically the Goldstone mode is responsible for long-range interactions between 
vortices in the superfluid phase, and at non-zero T for propagating wave excitations of local 
energy density known as second sound. 

In order to study the spectrum in the spin— | sector we need the Gor'kov propagator 
Q = A~ x . It can be written as 



A N 

-t^xy ly xy 

AT A 

ly xy -t^xy 



(16) 



where each entry is a 2 x 2 matrix in isospace. The notation indicates both normal (q(x)q(y)) 
and anomalous (g(x)q(y)) components: on a finite system A and A vanish in the limit 
j,j — > 0, and we recover the usual fermion and anti-fermion propagators as N and N 
respectively. The number of independent components is constrained by the identities N 22 = 
^21 = _ N* 2 , ^22 = —A^ and A 2 i = A\ 2) so that Q can be reconstructed using just two 
conjugate gradient inversions of A. In addition, isospin and time-reversal symmetries dictate 
that after averaging over configurations the only independent non- vanishing components are 
ReiVn and Iul4 12 , which henceforth we will refer to simply as N and A respectively. 



C. Phenomenological Parameter Choice 

As previously discussed, in a number of dimensions d > 4, the continuum NJL model is 
non-renormalisable, which implies that the lattice model becomes trivial in the continuum 



limit 
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means one must choose a fixed lattice spacing a, corresponding to a fixed 

2 

inverse coupling constant (3 = ^2, by fitting the model's parameters to low energy, vacuum 
phenomenology. Employing methods outlined in j3| , we calculate the ratio between the pion 
decay rate /„. and the constituent quark mass, i.e. the dynamical mass of the quark in the 
chirally broken phase. By fitting to phenomenological values we may extract (3 as a function 
of the model's only other free parameter, the current quark mass m. Finally, calculating and 
fitting the pion mass allows us to fix m, and hence (3. We take advantage of the fact that 
a perturbative expansion in 1/N C is possible in four- Fermi theories by calculating quantities 
analytically to leading order in 1/N C , the Hartree approximation. Feynman diagrams are 
evaluated using staggered quark propagators defined on an x L t /2 A Euclidean blocked 
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lattice with periodic boundary conditions in spatial dimensions and anti-periodic boundary 
conditions in the temporal dimension. The integrals over loop momenta are evaluated as 
sums over momentum modes. 

Let us first calculate the gap equation, the fermion self-interaction, to leading order 
in 1/N C . For sufficiently strong coupling g 2 > g 2 the scalar auxiliary field a develops a 
spontaneous vacuum expectation value S, which in the chiral limit can be identified with 
the constituent fermion mass and is given by 

where p 2 = X^=o s * n2 P» * s ^ ne dimensionless squared loop momentum, N c and Nf are the 

number of flavors and colors respectively, and (3 = a 2 / g 2 is the dimensionless inverse coupling 

constant. In this section we shall distinguish the expectation value of the scalar field £ from 

the constituent mass, which to leading order in (1/N C ) we define as X* = X + m; elsewhere 

in this paper the two shall be used interchangeably and both denoted by X. 

/tt is calculated from the vacuum to one-pion axial-vector matrix element 

(p+ k/2)a- 1 

(o \JisM\iCj) = h^i5Y ^^^9^^> ka -+ ( 18 ) 

(p - k/2)a- 1 

Translating this diagram, noting that the pion-quark-quark coupling g nqq ~ g/^N c , and 
taking the k — > limit we find 

d 4 p cos 2p u 

(19) 



2N c Nt 

U " ( 27T ) 4 [P 2 + (S*a 




d 4 p 1 



2 



'-f ( 27r ) [p 2 + (S*a) 2 ]' 

In order to check that this form is sensible we choose to examine the continuum limit 
by extracting the leading order behaviour of (JT§j) as the dimensionless quark mass £*a is 
reduced to zero. This is done by the introduction of a dimensionless hyper-spherical cut-off 
S which splits the loop momenta into two regions, one with \p\ > 5 and one with \p\ < 5. 
As 8 — > 0, the inner region picks up the continuum behaviour, whilst the outer region yields 



the finite terms relevant to lattice perturbation theory 



3- 1 



gnoring these terms and taking 
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£*a — > we pick up a leading contribution f^/T,* ~ ^jN c Nf ln(<5/£*a)/47r 2 . Although 
this quantity is logarithmically divergent, we know that the integral between ir/2 and — 7r/2 
is both finite and independent of 5, and that the transition between the two regions of 
integration is smooth. This means that the In 8 term must be cancelled out by a similar 
term in the outer region and the continuum behaviour of /tt/S* is thus given by 



S* y (2tt) 2 X*a V ; 

This is the same behaviour found as S*/A — ► for the regularisation schemes employed 
in namely 3<i-momentum, 4<i-momentum, and real-time cut-offs as well as the Pauli- 
Villars scheme. 

By calculating (|19|) in the infinite volume limit, fitting f n to its experimental value of 
93MeV and £* to a reasonable 400MeV we are able to extract the dimensionless quark mass 
£*a = 0.557, such that to leading order in l/N c the lattice spacing a = ^OMeV)" 1 ~ 0.3fm. 
Solving the gap equation (JT7J) with this value for the mass we find that /3£a = 0.273. Using 
the identity S* = m + E we deduce a relationship between the bare quark mass and the 
inverse coupling: 

0^73 

A (0.557 -ma) K ' 

Finally, in order to fix the bare quark mass and hence the coupling, we need to fit one more 

phenomenological observable. Again following l^ . we calculate the mass of the pion m n by 

solving the self-consistent equation 

4- = 2N c N f ^Ll, (22) 

where I represents the integral 

d A p 1 



-| (2vr) 4 \f+ + (S*a) 2 ] [p 2 _ + (S*a) 2 ] ' (23) 
and p± = Ylt=o s i n2 (p=t^ m 7ro);y. Setting to a phenomenologically reasonable 138MeV and 
demanding that (}2T]) is satisfied fixes the bare mass to ma = 0.006 and the inverse coupling 
to (3 = 0.495. Table contains a summary of the fits made and parameters extracted. 



D. The Simulation 



We chose to perform the numerical simulation of the path integral defined by (J7J) in the 
partially quenched limit j = j = 0, for two reasons. The first is technical; in this limit 
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Phenomenological 
Observables Fitted 


Lattice Parameters 
Extracted 


S* = 400MeV 
U = 93MeV 
= 138MeV 


ma = 0.006 
(3 = 0.495 
a' 1 = 720MeV 



TABLE I: Summary of large- N c parameter fits. 



Pf 2 (,4) = detMM\ so that there is no obstruction to using the hybrid Monte Carlo (HMC) 
algorithm 2fJ, which is "exact" in the sense that there is no systematic error due to non-zero 
time step, a reassuring feature whenever a previously unexplored phase is to be studied. The 
second is pragmatic; it enables many values of j to be studied with a sin gle run at fixed g 2 
and fi, thus saving much computer time. Our previous studies in 2+ld jrsl \v\ used both 
partially quenched and full pfaffian simulations, and found essentially no difference in the 
results. 

The simulations were run on a variety of lattice sizes and values of fi. In a typical 
simulation, we generated 500 equilibrated configurations separated by HMC trajectories of 
mean length 1.0. The time step required to achieve an acceptance rate ~ 80% was found 
to decrease with increasing lattice size. For larger lattices, where the optimal time step 
was smaller than 1/25, we also tuned the number of colors N' c > N c = 8 used during 
the generation of molecular dynamic trajectories; in tuning the number of colors, which 
is "renormalised" by the discretisation of the trajectory, one can increase the acceptance 
without increasing the number of integration steps used. In the exact accept/reject step, 
and hence the integral over all configurations, N c = 8. Measurements were carried out on 
every second configuration, in all cases with j = j =^> j_ = 0. In the rest of the paper we 
will quote all results in terms of j + , henceforth referred to as j. 

The reader may be wondering why HMC simulations of the NJL model with \x ^ are 
possible, whereas they are well-known to be ineffective for QCD. In both cases the algorithm 
reproduces as path integral measure the positive definite detMM^. For QCD M describes a 
color triplet of quarks q, a color anti-triplet of conjugate quarks q c , thus permitting gauge 
singlet baryonic qq c bound states. The lightest such state can be shown to be degenerate with 
the Goldstone pion associated with chiral symmetry breaking, which means that baryonic 
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matter is induced into the ground state at an onset Uqhmc ~ 0(m n /2), in contradiction to 



the physical expectation [i qcd ~ 0(m nuc i eon /3) [2jj,|28|]. in NJL-like models, by contrast, in 
the large- N c limit the Goldstone pole is saturated by disconnected qq bubbles, which cannot 
contribute in the qq c channels [29(; hence the HMC onset occurs at the constituent scale 
Ho ~ £ as expected, resulting in a ground state of degenerate fermion degrees of freedom 
and hence a Fermi surface. For the continuum model studied here q and q c have identical 
quantum numbers and are hence indistinguishable; however the HMC approach also yields 
physically reasonable behaviour even in models where the chiral symmetry group is U(l) 
rather than SU(2)®SU(2) and this is not the case 2^ |. 

The observables defined as matrix traces in equations (|9llUll2j) may be estimated on a 
particular $ configuration by calculating the bilinear rfTA^r) using a conjugate gradient 
algorithm with a stochastic source vector 77 satisfying 77+77 oc 5 xy 5 pq ; during each measurement 
we used a set of 5 such vectors. The Gor'kov propagator Q is similarly calculated but 
using this time a point source. Special care is needed for susceptibilities, which contain 
contributions from both connected and disconnected fermion line diagrams, and hence must 
be calculated using both methods via expressions such as 

X = (^TA-^v^TA- 1 ^) - (^VArhf) 2 + £}<tr0 te r0£r>. (24) 

X 

We label these two contributions the disconnected and connected parts of x respectively. 



III. ZERO TEMPERATURE PHASE STRUCTURE 



A. Chiral Symmetry Restoration 

In this section we investigate the nature of the chiral restoration transition with our 
phenomenologically motivated parameter set, since the order of this transition is found to 
be sensitive to the choice of parameters employed Q|. We measure the chiral condensate 
(xx) an d baryon number density % on V = x L t lattices with L s = L t = 12, 16 and 
20, and for various chemical potentials [ia between 0.0 and 1.2. Henceforth all dimensionful 
quantities will be quoted in lattice units a = 1. The quantum average and statistical 
errors of the measured quantities are calculated using a jackknife estimate and the results 
extrapolated linearly to the infinite volume limit V' 1 — > 0. 
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Our results are presented in Fig. In order to compare the lattice data (points) with 
perturbative results, both (xx) an d are calculated to leading order in 1/N C (solid curves), 
corresponding to a mean-field theory in which the scalar field a = E on every dual lattice 
site and the auxiliary pseudoscalars 7Tj are exactly zero; in this limit the condensate is given 
by (xx) — = p' ( a )- At /i = 0, the large- N c solution predicts a non-zero condensate 
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FIG. 1: Chiral condensate and number density extrapolated to V^ 1 —> as functions of showing 
both the large- N c solution (solid curve) and lattice results (points). The dashed curve is propor- 
tional to the volume of the Brillouin zone bounded by the large- iV c Fermi surface. 

and zero baryon density, corresponding to a vacuum with broken chiral symmetry. As /i is 
increased the system passes into a phase in which chiral symmetry is approximately restored 
and matter begins to build up, with a pronounced "kink" at \i c ~ 0.75. In particular, ng is 
seen to be closely related to the volume of momentum space enclosed by the Fermi surface 
defined by the Fermi momentum kp given for free fermions by 

3 

sinh 2 fi = sin2 k Fi + (25) 

and plotted in the large-A^ c limit as the dashed line. This implies that in the continuum 
Ub oa k F as one would expect, whilst on a lattice % saturates at the value 1.0 as fi —>■ oo. 

The lattice data agree qualitatively with our analytic solution, although for these data 
both (xx) an d A*o are roughly 15% smaller, an effect we attribute to 0(1/ N c ) corrections. 
The transition at \i ~ fi c has the appearance of a crossover, and may thus be compatible 
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with a second order transition in the chiral limit m — * 0. If this is the case the region with 
riB > for « < fa could be associated with a "nuclear matter" phase, since if, as T — > 0, 
S falls below its vacuum value S for fj, < // c , the only possible physical agency is % > 
(a mixed phase of quark matter droplets in vacuum at constant pressure is not consistent 
with thermodynamic equilibrium unless there are at least two conserved quantum numbers 
present j^l). This behaviour is in marked contrast with that of the 2+ld model in which the 
transition is str ong ly first-order and baryonic matter has chiral symmetry restored at any 
density [17j, |l9|, |3l| . Finally, it is illustrative to convert these densities into physical units. 
The large- N c estimate a -1 = 720MeV translates the lattice point (/ia, n^a 3 ) ~ (0.65,0.072) 
into jx = 470MeV, n# = 3.5fm~ 3 . Bearing in mind that due to species doubling, in a cube 
of volume (2a) 3 \ describes two spin and four color components of a continuum fermion, we 
deduce a total physical density of 0.88iV c iV/fni -3 , to be compared with the nuclear matter 
onset fj, q ~ 320MeV, n q ~ 0.15A^ c fm -3 . 



B. Diquark Condensation 

The main purpose of this study is to determine the nature of the high density phase 
in which chiral symmetry is approximately restored. In particular, in order to explore the 
possibility of a U(l)#-violatmg BCS phase, we study the diquark order parameters ()12|) and 
their susceptibilities (114)) as functions of chemical potential. On their own, the susceptibilities 
are of limited importance. In conjunction with the Ward identity (fTK|) however, their ratio 

R=-^± (26) 

X- 

provides an important tool to distinguish between phases in which U(1)b symmetry is either 
manifest or broken. With manifest symmetry, and in the limit j —>■ 0, these two susceptibil- 
ities should be identical up to a sign factor, and the ratio should equal 1. If the symmetry 
is broken, however, the Ward identity predicts that the Goldstone mode X- should diverge 
as j — > and hence R should vanish. 

As stated in Sec. Ill Dl the disconnected terms of ()14j) are calculated via the use of multiple 
noise vector estimation, whilst the symmetry constraints discussed after (JTfij) imply that the 
connected terms are given by 

xT = \{- [\Gn\ 2 + l&il 2 + |£ 33 | 2 + |£ 43 | 2 ] ± [|£i3l 2 + Ifel 2 + l&il 2 + l&il 2 ]) , (27) 
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evaluated between a random point source and the point x. 

The susceptibilities are measured and R calculated on the aforementioned lattice sizes 
and for various values of //, with the diquark source j varying from 0.1 to 1.0 during each 
set of measurements. It is interesting to note that although in most cases the disconnected 
contributions are found to be consistent with zero, in the low fi phase with large j, yl* s 
can be up to 10-20% the magnitude of x+ n - 111 contrast with the NJL model in 2+ld [19], 
therefore, we cannot assume that \+ — X+ n - 

An interesting empirical observation is that whilst the observables measured in Sec. IIII Al 
were found to scale linearly with the inverse volume of the lattice, observables in the di- 
quark sector appear to scale linearly with the inverse temporal extent, corresponding to the 
temperature of the system. Accordingly, the ratio R is extrapolated linearly to the limit 
L^ 1 — ► 0, which as in is found to give a reasonable description of the data. An example 
of the quality of the fits is illustrated in Fig. EJ 
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FIG. 2: Susceptibility ratio R at /j, = 0.0 and 1.0 on various lattice sizes. 



Figure El shows this extrapolated data plotted as a function of j for various values of 
fi, as well as the results of a linear extrapolation to j — > 0. One immediately notices that 
whilst this extrapolation appears plausible for j > 0.3, the data at lower values of j in 
the high \x phase diverge rapidly from the linear trend. The fact that this effect increases 
systematically with increasing \i and decreasing j suggests that its origin is some systematic 
effect not considered thusfar; indeed, the study presented in Sec. EU shows this to be due 
to residual finite temperature effects. For this reason, we believe that we are justified in 
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FIG. 3: R extrapolated to j — > for various values of (j,. 



disregarding data with j < 0.2 when taking the j — > limit. With this omission Fig.Elshows 
that for fi = a linear fit is consistent with a ratio of R ~ 1, corresponding to a manifest 
baryon number symmetry as one would expect in the vacuum. At /i = 1.0, however, R w 
suggesting that U(l)# symmetry is broken. 
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FIG. 4: (qq+) extrapolated to j — > for various values of /i. 



For more direct evidence of diquark condensation we measure the order parameter defined 
in (|12j) . Again, these data are extrapolated linearly to the limit L^ 1 — > with the quality 
of the fits being good. Figure |U shows the extrapolated values of (qq + ) plotted against j 



16 



for various values of fi. Fitting a quadratic curve through the data with j > 0.3, one can 
clearly see that the high /i, low j data again disagree with the curve; again ignoring these 
points, the data are extrapolated to j — > 0. For /i = we find no diquark condensation as 
one would expect, but as \i increases from zero, so does (qq+). Together, the observations 
that lim.j_ R = and linx^o (QQ+) 7^ support the existence of a BCS superfluid phase at 
high chemical potential. 
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FIG. 5: A comparison between the diquark and chiral order parameters as functions of [i. 



Finally, (qq+) is plotted as a function of \i in Fig. along with the previously presented 
result for (xx)- Although there is clearly a transition from a phase with no diquark con- 
densation to one in which the diquark condensate has a magnitude similar to that of the 
vacuum chiral condensate, this transition is far less pronounced than in the chiral case. 
(qq+) increases approximately as /z 2 , but eventually saturates as fi approaches 1.0 and even 
decreases past fi ~ 1.1. This behaviour is directly related to the geometry of the Fermi 
surface for a system defined on a cubic lattice, the area of which we have calculated in the 
large- A^ c free fermion limit from and have plotted as the solid curve. The method of this 
calculation is sketched in Appendix El In the continuum, (qq+) should continue to rise as 
[i 2 , such that the curvature d 2 (qq+) /dfi 2 is positive, in contrast to the behaviour observed 
in simulations of two color QCD, in which there is no Fermi surface and U(1)b breaking 



proceeds via Bose-Einstein condensation 



32 



33 



leading to d 2 (qq+) /dfi 2 < 0. 



The apparent weakness of the transition at intermediate \x is related to the fact that at 
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these chemical potentials, the value of R\j_> interpolates between the two extremes of 
and 1. This is counterintuitive, since it suggests a partially broken symmetry, even at j = 0. 
It may be that this is a side-effect of the chiral transition being a crossover, since there is 
no sharp point at which a large Fermi surface is created. It is also possible, of course, that 
this behaviour for intermediate \x is merely an artifact of our poor control over the j — > 
extrapolation. 



IV. THE QUASIPARTICLE DISPERSION RELATION 
A. Spectroscopy in the Fermionic Sector 

In this section we study the dynamics of the model's fermionic excitations, which as 
in the original BCS theory can be viewed as quasiparticles with energy E relative to 
the system's Fermi energy Ep. In a traditional Fermi liquid, these can be identified with 
particle excitations above the Fermi surface, and hole excitations below, both of which can 
have energies arbitrarily close to E — 0. In a superfluid system, however, the particles and 
holes mix, and energies of the lowest-lying excitations are separated from zero by a BCS 
gap A, which in analogy to the chiral mass gap in the vacuum Eq, can be viewed as an 
effective order parameter for the system. One advantage of this parameter is that unlike the 
diquark condensate, A can be directly related to a macroscopic thermodynamic property of 
the system, the critical temperature T c (3^ . In principle, being a spectral quantity it is also 
measurable in a color superconducting phase in QCD, where according to Elitzur's theorem 
one cannot define a local order parameter in a gauge invariant way [36]. 

The propagation of the quasiparticles is described by the Gor'kov propagator, defined 
in (J16|) . such that by analysing its momentum dependence we can map out a quasiparticle 
dispersion relation E(k), i.e. the energy spectrum as a function of momentum, and hence 
measure A. In particular, we have measured the time-slice form of both "normal" and 
"anomalous" propagators 

N(k; t) = Re (0, 0; x, t) e - ^*; (28) 

X 

A(k;t) = J2^ [A12 (0, 0; x, t)] e~ i%£ , (29) 

X 

on L x x l? y z x L t lattices with L x = 96, L y z = 12. This choice of having one spatial direction 
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much longer than the others gives the system a large number of modes with which to sample 
E(k), whilst minimising the extra computational expense of running on a larger volume. In 
particular, by setting k = (/c, 0,0) with k = 2nn/L x (n = 0, 1, 2, . . . , L s /4), the lattice 
fermions have 25 independent modes between and 7r/2 in the k x direction. Since the study 
presented in Sec. |V| shows that the diquark observables display little spatial dependence, 
there should be no detrimental effects from working with L x ^> L V}Z . 

Simulations were performed with L t = 16 and 20 at various chemical potentials using 
the same values of the diquark sources as in the previous section. L t = 12 data were also 
generated, but these prove to have too few time-slices over which to reliably fit the prop- 
agator, and are therefore neglected in our analysis. Again, approximately 500 equilibrated 
trajectories were generated per run, with measurement taking place on every other configu- 
ration. Two additional simulations were performed at fi = 0.0 and 0.8 with L t = 24, which 
after equilibration took approximately 5| and 16 CPU days respectively to generate 400 
trajectories on a 2.0GHz Intel Xeon processor. 




t t 

FIG. 6: Normal and anomalous propagators measured on a 96 x 12 2 x 16 lattice in both the chirally 
broken and restored phases. 

Some example propagators in both the chirally broken and restored phases are plotted 
in Fig. El At /i = 0.0 the normal propagator is non-zero for all t, whilst at fi = 0.8 
it approximates zero on even time-slices, which reflects the fact that with manifest chiral 
symmetry the N QO and N ee components of the standard staggered fermion propagator vanish. 
The anomalous propagator is zero on all odd time-slices for all values of fi. 
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To map out the dispersion relation for each value of /i, the energy is extracted by fitting 
the propagators to 

N(k,t) = Ae~ Et + Be- E(Lt -V if t = odd 
N(k, t) = if t = even 

and 

A(k,t) = C(e~ Et - e - E{Lt ~^) if t = even 
A(k,t) =0 if t = odd, 

where A, B and C are kept as free parameters, as is the energy E. Although, as expected, 
the values of E extracted from the two propagators are found to be consistent, we choose to 
use those extracted from (J31|) to map out our dispersion relation, since having one less free 
parameter than (|30p. the fits to this form are found to be of a higher quality. Some example 
fits are illustrated in Fig. [7| 



(30) 



(31) 



"i 1 1 1 1 1 r 



A = 0.326(e-°- 547 * - e -o.647(z«-t)) 
t e [4, 12] r7 dof ~ 0.33 




10 12 14 16 



FIG. 7: Normal and anomalous propagators with j = 0.3 and k = tt/4 with their fitted curves 
measured on a 96 x 12 2 x 16 lattice at n = 0.8. 



Figure |H1 shows plots of the free parameters in (jSOjl and (pTTj) at [A — 0.8, extrapolated to 
T — > 0; in turn these are extrapolated to j — *■ 0. Quadratic polynomial curves are fitted to 
the coefficients A(k), B(k) and C(k), whilst the energy E(k) is fitted linearly. As with the 
extrapolations of (qq+) and R, these appear to smoothly fit the data except for those with 
low values of the diquark source j, for which the discrepancy discussed in Sec. IIIIBI persists. 
Again, points with j < 0.3 are ignored for the purpose of the extrapolations; we rely on the 
conclusions of Sees. |V|and|VT]to justify the omission of these data from our analysis. 
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FIG. 8: Zero temperature propagator parameters at fj, = 0.8 and various values of k. The solid 
curves show extrapolation to j — ► 0. 

B. The Vacuum Dispersion Relation 

Before considering a lattice system with a Fermi surface, we review the nature of the 
dispersion relation in the familiar case of the vacuum. With fi = 0, time reversal is a 
good symmetry of the lattice and the coefficients A and B become identical such that ()30|) 
reduces to its usual form \N(k, t)\ = A(e~ Et + e~ E ^ Lt ~ t ^'). This can be understood physically 
by noting that the vacuum spectrum appears identical to both particles and antiparticles, 
and hence to both the forward- and backward-moving parts of N. In agreement with this, A 
and B are found to be equal, within errors, for all three values of L t . Figure El illustrates the 
energies extracted from A(k, t), extrapolated to zero temperature through L t = 16, 20 anc 
24 and then to j — > 0, which results in the familiar lattice dispersion relation E(k, fi = 0) 

3 

sinh 2 E = a 2 sin 2 ki + Sq, 
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(32) 



21 



k 



FIG. 9: The vacuum dispersion relation E(k) extrapolated to T —* and then j — > 0. 

where a is a constant which allows for the renormalisation of the speed of light by thermal 
effects; in practice its fitted value is close to one as expected. The energy gap at k = can be 
identified with the vacuum fermion mass, from which we learn that S = 0.351(6). As k is 
increased, the dispersion relation is approximately quadratic for small k/T, (as expected for 
a non-relativistic particle), until discretisation effects dominate its form and the periodicity 
of (p?^ causes it to level off as k — > n/2. 

C. Measurement of the Gap 

In this section we study the dispersion relation at /i = 0.8, with the aim of observing 
the BCS gap A. First, however, we study the momentum dependence of the other free 
parameters fitted from the forms fK~>OI) and Ijfllj) . Figure ITUl illustrates the values of A, B and 
C extrapolated first to T — ► and then to j — ► 0, and plotted as functions of momentum k. 

The coefficients A(k) and B(k) represent the amplitudes of forward- and backward- 
moving spin- 1 propagation, which due to our choice to study the antiparticle propagator 
iV = Nu = Gsi ~ (q(x)q(y)), correspond to hole and particle excitations respectively in the 
limit that j — > 0. For small momenta, corresponding to excitations deep within the Fermi 
sea, propagation is dominated by hole degrees of freedom. As the Fermi momentum is ap- 
proached, particles are easier to excite and become the dominant contribution as k — >• f . 
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FIG. 10: The propagator coefficients A, B and C for fj, = 0.8 extrapolated to T — > and then 

To degrees of freedom with = as in the vacuum the background appears the same 
to both particles and holes such that A(kp) = B{kp). For this reason, in analogy with 
the coefficients of filled and unfilled states in the original BCS theory |l6[, the intercept of 
the curves of A(k) and B(k) defines the Fermi momentum for the interacting theory in the 
T — > limit. In the anomalous sector, the coefficient C(k) is approximately zero deep within 
the Fermi sea, but becomes non-zero (even in the limit j — > 0) in a broad peak about the 
Fermi momentum, which is a sign of particle-hole mixing in this region. This is in contrast 
with similar measurements in NJL 2 +i jl^t - 

The left-hand panel of Fig. ^2 illustrates the lattice dispersion relation E{k) at [l = 0.8 
extracted from (}3*T|) and extrapolated first to T — > and then j —>■ (points), compared 
with that of free massless fermions on the lattice, parametrised by 

E(k) = -/i + sinh _1 (sinA;). (33) 

For this free dispersion relation, there are two distinct branches, one for hole excitations 
below the Fermi surface for which E reduces with k, and one for particle excitations for 
which E increases with k. We should note at this point that the crossover between these 
regimes, kp = sin _1 (sinh0.8) ~ 0.3487T, is consistent with the intercept of A{k) and B{k) in 
Fig. El to within the precision allowed by the momentum resolution. By contrast the lattice 
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FIG. 11: The lattice dispersion relation and typical free fermion dispersion relation at /i = 0.8. In 
the right-hand panel the hole branch is plotted as negative. 

data display no evidence for two distinct branches, which is another signal of particle-hole 
mixing. More importantly, at no point do the lattice data pass through zero, but between 
7r/3 < k < 177r/48 (again consistent with the free-field hp), they have a minimum of 
E(k) = 0.053(6); this is the BCS gap A. Comparing this to our measurement of the vacuum 
fermion mass in Sec. II V Bl we find the ratio 

A(/i = 0.8) 



0.15(2), 



(34) 



which assuming a fermion mass of 



with the analytic predictions of 
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OOMeV implies that A(/i = 0.8) ~ 60MeV, consistent 



This may be viewed more graphically in the right-hand panel of Fig. ^2 where excitations 
below Ep have been plotted as negative. This makes the free fermion dispersion relation 
a smooth curve that passes continuously through E = at the Fermi momentum, and is 
similar in nature to that of lattice four- Fermi models with no BCS gap [rol . [^J. For our 
data, however, this plot introduces a discontinuity of 2A which looks exactly like that of a 
traditional BCS superfluid in the continuum. 



D. fx Dependence of the Gap 

Having systematically investigated the dispersion relation and extracted the gap at one 
value of chemical potential, it would be illuminating to repeat the analysis of the previous 
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section for a range of chemical potentials in the BCS phase. However, generating data with 
L t = 24 in the chirally restored phase is a CPU intensive task, taking O(20) CPU days on a 
fast desktop PC for each value of //. The reason this is so much more expensive than in the 
chirally broken phase is that the rate of convergence of the conjugate gradient subroutine 
used to invert M^M in the generation of our {$} configurations is related to the magnitude 
of the diagonal components of M, which are in turn proportional to the constituent quark 
mass > fi c ) ~ 0. 

For this reason we utilise the data generated with L t = 16 and 20 and approximate the 
T — > limit by extrapolating through these data and assigning a conservative estimate 
for the error; the j — > extrapolations may then be carried out as before. Although of 
little statistical significance, it should be noted that the dispersion relation at fi — 0.8 
extracted using this method is consistent, for all values of k, with that extracted using the 
full statistical treatment of the previous section. 
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FIG. 12: Lattice dispersion relations at various values of fi. 

Figure El illustrates a selection of these dispersion relations for various chemical poten- 
tials. An interesting point to note is that unlike the diquark condensate, there is an upper 
bound of /i ~ sinh~ 1 (1.0) ~ 0.88, above which we cannot extract an estimate of the gap 
using the method outlined in Sec. IIV Al This can be understood by considering the nature of 
the lattice Fermi surface for free fermions in the infinite volume limit, parametrised by ()25|) 
and plotted at two values of fi in the large- N c limit in Fig. El For kpi <C vr/2, corresponding 
to small values of the small-angle identities imply that the surface is approximately spher- 
ical and the momenta sampled herein, emphasised on the k x axis in Fig. El are sufficient to 
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FIG. 13: The zero temperature, infinite volume Fermi surface at two values of chemical potential. 



probe either side of the surface and detect the presence of a gap. For /i > sinh - (sin vr/2), 
however, the reduced rotational invariance of the lattice dominates, and (J25|) has no real 
solution with k = (0 < k x < |, 0, 0). It is for this reason that the curve at /i = 0.9 in Fig. 
shows no minimum and we cannot extract a value of A. Whilst there is almost certainly a 
gap present, as the maximum of (qq+) lies between 1.0 < // < 1.1, to detect its presence via 
the dispersion relation one is required to sample momenta along e.g. the more complicated 
diagonal path k = (k, k, k) with < k < tt/2, illustrated in the right-hand panel of Fig. 
Because of the large number of spatial modes required to sample this path with sufficient 
resolution, such a study is computationally beyond our current capability. 
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FIG. 14: Estimates for the gap A compared with {qq+) for a range of jJL. 



Finally, A is plotted as a function of /i in Fig. El and compared with the value of the 
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diquark condensate. Although the error bars on our estimated values are fairly large, we see 
little evidence of any \i dependence once the gap becomes non-zero. In fact, a least-squares 
fit to A = constant, denoted by the horizontal bar, has a chi-squared value of only 0.33 per 
degree of freedom. In combination with Fig. Fig. provides qualitative support for a 
simple-minded picture in which only diquark pairs within a shell encasing the Fermi surface 
of thickness 2A, independent of /i, participate in the pairing, resulting in a condensate 
(qq+) oc A/z 2 . 



V. FINITE SPATIAL VOLUME EFFECTS 



The conclusions of Sees. IIHI and HV] that the high-// phase is one with both (qq+) ^ 
and A ^ 0, both rely on the discarding of data with j < 0.3, since results in the diquark 
sector with these small diquark sources disagree with the trends at higher j. In order to 
be able to trust our j — > extrapolations it is necessary, therefore, to justify this omission, 
especially since it is the data in the region of this limit that have been discarded. 

We have previously argued that the discrepancy at low-j could be due to finite size 



effects 22|, since whilst our simulations were performed on lattices with L s < 20, variational 
studies of the Nf = 2, N c = 3 continuum NJL model at zero temperature in a finite spatial 



mj~ 25 lattice spacings) is 
We have argued further 



volume show that with no diquark source, a spatial extent of 7: 
required before the model approximates its infinite volume limit 
that the source of these finite size effects is due not to the realisation of an exact Goldstone- 
mode, but to the difficulty of representing a thin shell of states about the Fermi surface 
which contribute to diquark condensation on a discrete momentum lattice. Whilst this is 

n 

supported in part by the finite size study presented in [22j, in which (qq+) displays a non- 
monotonic dependence on L s , it should be noted that the magnitudes of these fluctuations 
are less than 1% for L s > 8 and all values of j, much smaller than the approximate 30% 
suppression of (qq+) at j = 0.1 seen in Fig. HJ In fact, the diquark condensate shows little 
notable L s dependence even prior to extrapolation to T — > 0. In Fig. HHJ (qq+) is plotted 
as a function of j at fi = 0.8 with L t = 12 and various spatial extents, including the results 
of a large simulation with L s = 30. This shows that the diquark condensate displays no 
significant size dependence at any j. The same lack of any significant spatial size dependence 
is found for both (qq+) with L t = 16 and 20 and for E(k) with L t = 16 and 20 within the 
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FIG. 15: Diquark condensate at \i = 0.8 on L x x Ly z x Lt lattices with Lt = 12 and various L 



allowed momentum resolution. The fact that these observables remain unchanged at low 
j, even when changing the spatial volume by up to a factor of 2.5 3 ~ 16, makes it hard to 
believe our previous suggestion that the suppression at j < 0.3 is due to finite size effects. 
Instead, we propose an alternative suggestion in the following section. 



VI. NON-ZERO TEMPERATURE 

In previous reports of this work 39, we have suggested that an obvious extension 
would be to carry out simulations at non-zero temperature, with the aim of measuring the 
critical temperature of the superfluid phase, T c . The non-relativistic BCS theory predicts 
the relation between this parameter and the magnitude of the gap to be [l(| 

£ = 1-764. (35) 

Furthermore, it has been shown that this relation holds for relativistic color superconduct- 
ing systems in weakly coupled QCD with two flavors j^]. Although such weak coupling 
predictions may be trusted only at asymptotically high densities, naively applying (|35j) to 
our measurement of A(fi = 0.8) = 0.053(6) suggests that T c ~ 0.03a -1 , such that at this 
chemical potential and in the limit j — > 0, one should observe a superfluid phase only when 
the temporal extent is greater than about 35 lattice spacings. The fact that we observe a 
BCS phase, even though our simulations were performed on lattices with temporal extents 
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much smaller than this relies on our performing measurements with j ^ and then extrap- 
olating in the correct manner. Setting j ^ has the effect of making condensation more 
favourable, which suggests that at fixed j there could be a crossover at some pseudo-critical 
temperature, T c (j), separating a region where diquark condensation is suppressed by ther- 
mal fluctuations and one in which it is not. One would expect the effect of increasing the 
source would be to increase T c (j) such that diquark condensation can be observed on lattices 
with smaller temporal extents. If one can successfully extrapolate to zero temperature first, 
a j — > extrapolation should then be possible. 

This causes a problem if one wishes to determine the value of the condensate at a partic- 
ular chemical potential and non-zero temperature, since it is possible that at temperatures 
close to T c , T c (j) crosses the temperature of the lattice over the range of j studied. Figure Hoi 
illustrates the diquark condensate measured at fi — 0.8 on lattices with various temporal 
extents corresponding to various non-zero temperatures, as well as the "T = 0" curves for 
fi = 0.0 and 0.8 from 12 4 , 16 4 and 20 4 lattices, as plotted in Fig. HJ The L t = 4 results lie 
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FIG. 16: Diquark condensate at fi = 0.8 on various lattices at non-zero temperature. 

well below the T = curve for all values of j, suggesting that the temperature is high enough 
to suppress condensation for the entire range of j. A quadratic extrapolation through these 
points is very similar to the zero temperature vacuum fit from Fig. H] and is consistent with 
(qq)\j_+ = 0. As the temperature of the lattice is decreased, however, the data at higher 
j are no longer suppressed and an extrapolation through all j is no longer possible. This 
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can be seen more clearly in Fig. El where we focus on the data measured on a 30 3 x 12 
lattice. Whilst an attempted fit through all values of j appears of poor quality, by choosing 
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FIG. 17: Diquark condensate at /i = 0.8 an a 30 3 x 12 lattice. 

a suitable point to separate the low- and high-j data one can fit the two regions successfully. 

Whilst this implies that establishing the critical temperature of the superfluid phase is 
not as simple as we initially believed, it does provide an alternative explanation of the 
suppression in (qq+) at T = for j < 0.3. The fact that the curve fitted through the 
L t = 12 data over j G [0.6 : 1.0] agrees with that of data extrapolated to T — > and fitted 
over j G [0.3, 1.0] suggests that these curves represent not the correct infinite volume limit, 
as previously suggested, but the correct zero temperature limit of the j — > extrapolation. 
Whilst a linear T — > extrapolation is sufficient to reach this curve for 0.3 < j < 0.5, the 
data at j < 0.3 must be suppressed too much for such an extrapolation to be sufficient, i.e. 
the temperature of the lattice must be too high compared with T c (j < 0.3). This justifies in 
retrospect the discarding of the low j data in the extrapolations of Sees. IIHI and IIVI which 
form the basis for the claimed superfluid ground state. 

VII. NON-ZERO ISOSPIN CHEMICAL POTENTIAL 

In previous sections, we have demonstrated numerically that the ground-state of our 
two flavor model at non-zero /i and low T is that of a BCS gapped lattice superfluid. 
The pairing mechanism in this system, between quarks of opposite momenta and isospin, 
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and thus analogous to 2SC condensation between u and d quarks in QCD, is particularly 
energetically favourable. If the system were non-interacting, it would cost no energy to 
create a pair of quarks at the common Fermi surface of the degenerate flavors, such that 
when the attractive interaction is restored, the system is expected always to be unstable to 



4l|. 



diquark condensation 

In nature, however, the Fermi momenta k F and k F for up and down quarks respectively 
are expected to differ. A simplistic argument outlined in [42], simplified still further here to 
describe a two flavor system, suggests that in compact stellar matter k F should be less than 
k F . For massless non- interacting matter with baryon chemical potential /jlb = 400MeV, an 
electron chemical potential /i e = 89MeV is required to enforce both charge neutrality and 
chemical equilibrium under weak interactions. Together, these two conditions determine all 
the chemical potentials and Fermi momenta: 

kp = \i u = \ib — §£*e — 355.5MeV, 

kp = ^ d = fx B + Ifie = 444.5MeV, (36) 
k*, = n e = 89MeV. 

Here we use the term baryon chemical potential in the context of the NJL model, where 
baryons are identified with quarks and \ib = §(/■*« + H-d)- The effect of separating the free- 
particle Fermi surfaces of pairing quarks should be to make the superconducting phase less 
energetically favourable, and should prove a good method to investigate the stability of the 
superfluid phase. 

Such a study was applied to the 2SC color superconducting phase in a mean field four- 
Fermi model in j43j. Similar to results for superconductors in the presence of a magnetic 
field 0, 45], when the free field Fermi surfaces are separated by only a small amount the 
ground-state of interacting matter remains superconducting with degenerate Fermi surfaces 
for the pairing partners. At some critical free fermion separation, however, the system is 
found to go through a first order transition to a gapless Fermi liquid with two separate 
surfaces. Unlike normal superconductors, however, the size of the gap increases slightly 
under small flavor asymmetries, an effect attributed to the model's color structure extracted 
from one gluon exchange in QCD. This analysis has also been extended further to include 
systems in which the Fermi surfaces k F and k d F are separated not only by an isospin chemical 



potential, but a fixed momentum q 



41 1 . In such a system, when \k F — k F \ and q are 



sufficiently large that the Fermi surfaces cross, Pauli blocking implies that 2SC becomes 
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unstable with respect to a state in which diquark condensation occurs only at a ring of 
states close to the intercept of the surfaces. The state has both broken translational and 
rotational invariance in which the diquark pairs have non-zero total momentum; in analogy 
with such phases in electron superconductivity this is known as the LOFF phase Q, |^|. 

In the lattice NJL model, the pairing quarks of opposite isospin are represented by the 
two components of the staggered fermion field 



/ Xi(aO ' 


H 


' u{x) \ 


\ X2O) I 




l d(x) ) 



^ \X2(X)) \d(x)) , ^ 

hereon referred to as "up" and "down" quark flavors. The Fermi surfaces of the pairing 
partners can be separated by directly allocating them different chemical potentials, fi u and 
fid, equivalent to having simultaneously non-zero baryon chemical potential fiB = \ + Ha) 
and isospin chemical potential fij = | (fi u — fid)- Although this definition implies that fi u > 
fid, which is contrary to the conclusions of the argument outlined above, this notation has 
been chosen to be consistent with the analytic studies of [48] and 3] (since the NJL model 
does not include weak interactions and is therefore isospin invariant, the labels u and d are 
interchangeable). In the physical context of compact stars, the two scales should be ordered 
fiB 3> Hi, since the simple argument leading to (jHEJ) predicts that | \fi u — fid\ ~ 0.1/!^. 
With this introduction, the fermion kinetic operator M, defined previously in (JHl) becomes 

- (38) 



i€{x)n(x).T pq ) 



(x,x) 

Unfortunately, this means that the proof that det M is real and positive presented in |Tj| is 
no longer valid, which can be seen be noting that e.g. T 2 (e T3flI )T2 = e~ T3tl1 7^ (e T3M/ )*, and the 
identity T2MT2 = M* no longer holds. Although this would not cause the simulation to fail, 
since we use det M^M as our fermionic measure, the fact that this choice is the sole reason 
the algorithm would work implies that non-trivial interactions between \ an d C quarks will 
be introduced which could cause the argument of [29] to break down. Instead, we choose to 
simulate in the quenched isospin limit in which det M^M is calculated with fij — in the 
HMC update of the bosonic fields, whilst (|3~8"j) is used in A during the measurement of the 
observables. 
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Before we discuss setting /// 7^ in the superfluid phase, let us examine the effect this 
introduction has on the chiral phase transition. An analytic study of the NJL model has 
shown that introducing a small /// causes the chiral phase transition to split into two, 
one transition for the condensate of each quark flavor j^. Figure illustrates lattice 
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FIG. 18: (uu) and (da) condensates for various hb and Hi on a 12 4 lattice. 



measurements of the up and down quark condensates 

„, W 4^ '/J 1± AaA (39) 

as functions of fiB for various \ij measured on a 12 4 lattice. Although these results are 
measured on only one volume, the speed of these simulations means that is is possible to 
gain fine resolution in yU#. Consistent with the predictions of the two transitions, which 
are coincident in the limit that ///—>• 0, separate as \ii is increased. This can be understood 
by noting that for fixed /ij the chemical potential of the up quark is larger than that of 
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the down, such that as //# increases, \i u reaches the critical chemical potential first. It is 
not clear, however, why the curve of the up condensate deviates from the \ii = solutions 
more than that of the down. This effect, not predicted in could be some finite volume 
artifact, or a result of the quenched approximation. As an aside, it has been argued that 
the observation of two transitions is an artifact of the diagonal flavor structure of the NJL 
model with broken chiral symmetry and would not be observed in nature. In particular, the 
introduction of an instanton-motivated flavor mixing vertex with even a weak coupling is 
shown to restore the single transition 

In the diquark sector, we relabel the order parameter (qq+) defined in (fT2"|) as (ud), to 
emphasise the fact that condensation occurs between quarks of different flavors. Figure HH] 

0.35 
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FIG. 19: (ud) condensate as a function of j for various /i/ with = 1.0. 



illustrates the (ud) condensate measured at \lb — 1.0 on 12 4 , 16 4 and 20 4 lattices as a 
function of j for various values of /ij. Results are extrapolated to T — > as before. As 
expected, the effect of significantly increasing fij for fixed j is to suppress condensation, an 
effect that is more pronounced at smaller values of j. Less straightforward to understand, 
however, is that at values of j above this suppression (ud(j)) appears to increase slightly 
with increasing /ij. By analogy with 43], this could be due to the non-trivial 0(a) color- 
mixing terms of the staggered quark action, but again could be due to the crudeness of the 
quenched approximation. 

As with our results at non-zero temperature, because the effect of increasing the diquark 
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source is to make the superfTuid more robust to setting /if 7^ 0, it is not yet clear how to 
determine the critical isospin chemical potential in the limit j — > 0. 

VIII. SUMMARY 

To our knowledge this is the first systematic non-perturbative study of a 3+ld relativistic 
system (albeit one with a Lorentz non-invariant cutoff) which has a Fermi surface. Our goal 
has been to study such a system with phenomenologically reasonable parameters, but with 
the main focus on characterising the high-density ground state. The principal result is that 
in the limit T —>■ the Fermi surface is unstable with respect to condensation of diquark 
pairs in the scalar isoscalar channel, leading to a ground state characterised by a U(1)b- 
violating order parameter (qq+). The resulting energy gap A which opens up at the Fermi 
surface is approximately 15% of the constituent quark mass scale, in good agreement with 
self-consistent estimates made with similar models in continuum approaches. It seems likely 
that the model is thus a superfluid described in terms of orthodox BCS theory. 

Detailed quantitative agreement with further aspects of the BCS theory, such as the 
prediction (J33j) for T c , is hard to verify because of practical difficulties in simultaneously 
probing both T — > and j — > limits. Finite volume effects in this system are complicated 
to unravel because of their separate dependencies on L s and L t . As described in Sees. Mand 
IVI| for the first time we have been able to simulate a sufficiently broad range of L s , L t to be 
able to demonstrate that the apparent suppression of the order parameter for small j is due 
to our distance from the L t — > limit. Note, though, that the ratio L t /L s cannot be made 
arbitrarily large without introducing artifacts due to discretisation of A;-space, as shown in 
Fig. 4 of |22j. With the bare parameters used in this study we estimate that volumes greater 
than 40 4 will be needed for quantitative studies of T c . 

Finally, we have made the first exploratory study of a system with non-zero chemical 
potentials for both baryon number and isospin, which as described in Sec. IVIIl is a necessary 
precondition for an accurate description of the conditions prevailing inside a neutron star. 
It is amusing that the sign problem comes back to bite, restricting us to considering non- 
zero isospin density in the quenched approximation only. Nonetheless, we have been able 
to observe the expected suppression of ud pairing at the Fermi surface. Bearing in mind 
that phenomenology demands /i/ <C /xg, it is not ruled out that future simulations could 
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successfully unquench isospin by either reweighting or analytic continuation in the small 
parameter //j///b, much as the conditions in heavy ion collisions can be accessed by lattice 
QCD simulations at /ig = analytically continued in /i^/T. 
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APPENDIX A: THE AREA OF THE FERMI SURFACE 

In Fig. |S] we plot the area of the Fermi surface in the infinite volume limit as a function 
of chemical potential \x. Here we outline the method used to produce this curve. 

In general, one can calculate the area of any surface S by integrating 0=1 over that 
surface 

Surface Area = / Ida, (Al) 
Js 

where da is an infinitesimal two-dimensional element of S. If the ^-component of that 
surface is single- valued, one may turn this into a two-dimensional integral in the x-y plane 

f ld(T - If TJr dydy > ( A2 ) 

Js JJAOxdy 

where A is the area projected onto z = and 5a is the area of a parallelogram tangential 
to S that projects onto a square of area SxSy. In the limit that S — > d, is given by the 
absolute gradient of z(x, y) 



\Vz(x,y)\ 



( dz dz dz\ 
\dx : dy } dz J 



' r\ \ 2 /a \ 2 

oz \ I oz\ 



aJ + UJ +1 (A3) 



and the equality in (jA2j) becomes exact. 

In calculating the area of the Fermi surface, we consider the section of momentum space 
with < k X) y )Z < | in which the surface is a single-valued function kF z (kF x ,kF y )- From 
(|25|). we find that the fc z -component of the Fermi surface for free fermions at fixed /i is given 
by 



kp z {kp xi kpy) = sin 1 y C — sin 2 kp x — sin 2 kp y (A4) 
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where the constant C = sinh 2 /^ — £ 2 and the value of is taken from the large- A*^ 

solution of the gap equation (fTTj) in the infinite volume limit. The absolute gradient of this 
function is 




, cos 2 fc Fa . sm k Fx + cos 2 sm fc F 

Vfcf^fc^fcj-y) = A /1+ - 2 , w-, 77—^ - 2, x (A5) 



and the surface area is given by 



.4 



\Vk Fz (k Fx ,k Fy )\dk Fy dk Fx , (A6) 



which we evaluate numerically. In setting the limits of integration, we have several cases, 
determined by the value of C(/i): 

(a) For C < 0, (jA4|) has no real solutions and the surface area is zero. From the definition 
of C one can see that this corresponds physically to the chemical potential being 
insufficiently large to allow the existence of a sea of particles with mass S. 

(b) For < C < 1, the approximately spherical surface intercepts the k x and k y axes at 
sin -1 \/C and we have one region of integration. 

(c) For 1 < C < 2, the discretisation of space-time dominates and the surface no longer 
intercepts the axes. This is the situation depicted in the right-hand panel of Fig. IT31 
(|A4|) no longer has any real solutions with sin 2 k x + sin 2 k y < 1 or > 2. We must 
separate our integral, therefore, into two regions. 

(d) For 2 < C < 3, the surface is again approximately spherical, but with the centre at 
(f i f > f )■ I n this final region, saturation effects are dominant as the surface area of 
the sphere decreases with increasing /z until at C = 3 (i.e. fi ~ sinh" 1 1-32) the 
area reduces to zero and the lattice is saturated with fermions. 

The limits of integration are listed in table |H] 

The only limitation of this method is that the anti-periodicity of the Fermi surface in 
the region considered implies that Vk Fz diverges at the boundaries k z = and k z — |. 
Analytically, this divergence is cancelled in the integral by the infinitesimal size of da; and 
dy. In a numerical solution, however both 5x and Sy are non-vanishing and the integral 
diverges. To overcome this effect we introduce a small "buffer" about these boundaries to 
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k x 
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< o 


Area = 
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: 1 


: sin" 1 VC 


: sin -1 \J C — sin 2 k x 


1 : 2 


: sin" 1 VC - 1 


sin -1 VC — 1 — sin 2 k x : | 


sin" 1 VC - 1 : f 


: sin -1 \J C — sin 2 k x 


2 : 3 


sin" 1 VC - 2 : f 


sin -1 a/C — 1 — sin 2 /c^ : f 



TABLE II: Limits of integration in calculating the area of the Fermi surface for different values of 
C = sinhV- S 2 . 

stop the inclusion of divergent terms. The size of this buffer is then reduced until its effect 
on the solution is negligible. The buffer used to produce the curve in Fig. 03 is 10- 9 ; once 
the curve is evaluated, it is multiplied by an arbitrary constant (1/45) to allow it to be 
compared directly with the measured value of (qq+). 

APPENDIX B: THE VOLUME OF THE FERMI SEA 

In Fig. ^ the volume of the Fermi sea in the infinite volume limit is plotted as a function 
of chemical potential. This calculation, which is simpler than that for the area of the Fermi 
surface, is done by integrating 0=1 numerically over the volume bounded by (|A4|) 

Volume = J J J ld^d^d^, (Bl) 

where again the limits are determined by the value of C. These are listed in table IIHI 

As the integrand, unity, is well behaved over all k x , k y and k z , there is no need to 
introduce a buffer into this calculation. Once the curve is evaluated, it is normalised such 
that Volume(C = 3) = 1 to allow direct comparison with the large- N c prediction for ns- 
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